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Abstract 

A special-purpose, semi-analytical solution method 
for determining the stress and deformation fields in a 
thin laminated-composite cylindrical shell with an 
elliptical cutout is presented. The analysis includes the 
effects of cutout size, shape, and orientation; non- 
uniform wall thickness; oval-cross-section eccentricity; 
and loading conditions. The loading conditions include 
uniform tension, uniform torsion, and pure bending. 
The analysis approach is based on the principle of 
stationary potential energy and uses Lagrange 
multipliers to relax the kinematic admissibility 
requirements on the displacement representations 
through the use of idealized elastic edge restraints. 
Specifying appropriate stiffness values for the elastic 
extensional and rotational edge restraints (springs) 
allows the imposition of the kinematic boundary 
conditions in an indirect manner, which enables the use 
of a broader set of functions for representing the 
displacement fields. Selected results of parametric 
studies are presented for several geometric parameters 
that demonstrate that analysis approach is a powerful 
means for developing design criteria for laminated- 
composite shells. 

Introduction 

Cutouts in cylindrical shell-type components are 
unavoidable in the construction of aerospace structures. 
This fact is significant because the structural failure of 
these components usually begins near the cutout 
because of high stress concentrations that initiate the 
formation of cracks. Hence, a cutout can trigger a local 
failure at a load level lower than the global failure load 
of a corresponding shell without a cutout. As a result, 
preliminary-design sizing of a cylindrical shell with a 
cutout is often based on the magnitude of the stress 
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concentrations near the cutout. Therefore, an accurate 
assessment of the stress concentrations in a given shell 
subjected to various types of loading and support 
conditions is essential to the development of safe and 
reliable designs. Moreover, validated special-purpose 
analysis tools that enable rapid parametric studies 
would be very valuable to structural designers and for 
the development of new design criteria and design 
concepts. 

Several analytical, numerical and experimental 
studies have been conducted during the past sixty years 
to determine stress distributions in cylindrical shells 
with a cutout and subjected to various types of 
loadings; such as, axial tension and compression, 
torsion, and internal and external pressure. Pioneering 
analytical work was conducted by Lurie 1,2 to investigate 
the effect of axial tension and internal pressure, and 
shell curvature, on the stress concentrations around a 
circular cutout in the 1940s. Many years later, 

analytical studies were presented by Lekkerkerker, 3 
Van Dyke, 4 Ashmarin, 5 Murthy et ah, 6 Guz et ah, 7 and 
Van Tooren et al. 8 that further investigated the effects 
of various factors on the stress concentrations around a 
cutout in a cylindrical shell. Similarly, experimental 
investigations have been conducted by Tennyson, 9 
Starnes, 10 Pierce and Chou, 11 Bull, 12 and Zirka and 
Chernopiskii, 13 and numerical studies have been 
conducted by Liang et ah, 14 and Shnerenko and 
Godzula. 15 In 1964 and 1972, respectively, Hicks 16 and 
Ebner and Jung 17 summarized the results obtained from 
several of these previous studies and provided extensive 
lists of references related to this problem. Most of 
these previous studies are for isotropic cylindrical shells 
with a circular cutout. Only a few of these studies, such 
as those presented by Pierce and Chou 11 and by Murthy 
et ah, 6 address the effects of cutout shape (elliptical 
cutouts) on the stress concentrations. 

Mitigation of high stress concentrations by 
tailoring shell-wall thickness, material orthotropy and 
anisotropy, and cutout reinforcement are also important 
considerations in the design of aerospace structures 
made of lightweight composite materials. Likewise, the 
potential for using shells of non-circular cross section 
are relevant to fuselage-like structures. However, only 
a few studies have considered these effects. For 
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example, the influence of wall-thickness variation on 
the stresses in axially loaded composite cylindrical 
shells, without a cutout, has been investigated by Li et 
al. 18 Although numerous analyses exist in the literature 
on the analysis of shells with circular cross sections, 
only a few include non-circular cross sections. 
Sheinman and Firer 19 provided an analytical 
investigation of stresses in laminated cylindrical shells 
with arbitrary non-circular cross sections. More 
recently, Hyer and Wolford 20 21 studied the effect of 
non-circular cross sections on damage initiation and 
progressive failure in composite cylinders by 
employing the finite element method. 

The objective of the present study is to present a 
special-purpose analysis for a laminated-composite 
cylindrical shell with an elliptical cutout that can be 
used to rapidly, and parametrically, investigate the 
effects of shell curvature; cutout size, shape, and 
orientation; and ply lay-up on stress-resultant 
concentrations near the cutout. The analysis is 
applicable to thin-walled cylindrical shells with non- 
uniform wall thickness, a non-circular (e.g., oval) cross 
section, and subjected to tension, torsion, and bending 
loads as illustrated in Fig. 1. 



Fig. 1 Geometry, coordinates systems, and applied edge 
tractions for an oval cylindrical shell with an 
elliptical cutout and non-uniform wall thickness. 


To accomplish this objective, an overview of the 
analysis is presented first. Next, the boundary value 
problem is defined along with the kinematics and 
stress-strain relations used in the analysis. Then, the 
derivation of the equations governing the response and 
numerical procedure are described. Finally, selected 
numerical results for oval and circular cylindrical shells 


with either circular or elliptical cutouts and subjected to 
either tension, torsion, or pure-bending loads are 
presented. 


Analysis Overview 

The analytical approach used herein permits the 
determination of the pointwise variation of 
displacement and stress components. It is based on the 
principle of stationary potential energy, but utilizes 
local and global functions that are not required to 
satisfy the kinematic boundary conditions directly. 
Thus, the choice of local and global functions is not 
limited by a particular type of kinematic boundary 
condition. The kinematic boundary conditions are 
imposed by employing the Lagrange multiplier method. 
Both local and global functions are used, in contrast to 
the traditional approach, to enhance the robustness of 
the analysis method. In particular, the local functions 
are used to capture rapidly varying stress and strain 
gradients and local deformations near a cutout. Toward 
that goal, Laurent series are used for the local functions 
and are expressed in terms of the mapping functions 
introduced by Lekhnitskii. 22 Fourier series are used for 
the global functions and are used to capture the overall 
deformation and stress fields. The kinematic 
admissibility requirements on the local and global 
functions are relaxed by defining that the edges of the 
shell are supported by extensional and rotational 
springs. Zero-valued displacement and rotation 
kinematic boundary conditions are enforced in an 
indirect manner by specifiying values for the spring 
stiffnesses that are large compared to the corresponding 
shell stiffnesses. This approach effectively yields a 
prescribed kinematic boundary condition in the limit as 
the relative stiffness of the spring becomes much 
greater than the corresponding shell stiffness. 
Similarly, values for the spring stiffnesses can be 
selected that correspond to a given uniform elastic 
restraint along an edge, similar to that provided by an 
end-ring. This capability is important, and useful, 
because in some test fixtures or actual structures the 
edge supports may not be stiff enough to simulate a 
fully clamped boundary condition or flexible enough to 
simulate a simply supported boundary condition. 

As suggested by Li et al. 18 and Sheinman and 
Firer, 19 nonuniform wall-thickness variations of a shell, 
which lead to non-uniform laminate stiffnesses, are 
represented by using trigonometric series. Specifically, 
nonuniform shell-wall thickness is represented in the 
present study by perturbing the ply thicknesses with a 
function that is periodic in either the longitudinal or the 
circumferential direction. The variation in wall 
thickness is accounted for by adjusting the lamina 
properties, resulting in nonhomogeneous in-plane and 
bending stiffness matrices. The nonuniform shell 
curvature associated with a noncircular cross section is 
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represented by using trigonometric series for the 
coordinates of an oval-cross-section shell reference 
surface. 23 The aspect ratio, or out-of-roundness, of the 
cross-section is represented in the analysis by using an 
eccentricity parameter. 

In the derivation of the equations governing the 
response, the total potential energy consists of the 
elastic strain energy of the shell, the elastic edge 
restraints and the potential energy of the applied loads. 
The conditions that may arise from the choice of 
displacement approximations without any kinematic 
restrictions are treated as constraint equations, and the 
potential energy arising from constraint reactions is 
invoked into the total potential energy through the use 
of Lagrange multipliers. The equations governing the 
shell response are obtained by enforcing the 
requirement that the first variation of the total potential 
energy vanish. The evaluation of the area integrals 
appearing in the potential energy are achieved 
numerically by using a basic quadrature method in 
conjunction with standard triangulation of the entire 
domain described by Shewchuk. 24 Solution to the 
equations governing the response are obtained by using 
a standared Gaussian elimination procedure, which 
yields the generalized displacement coefficients and, 
thus, the stress and strain fields. The accuracy of the 
analysis depends on the number of terms used for the 
functional representation of the displacement fields. As 
the number of terms increases, the results converge to 
the exact solution. 


Representation of Shell Geometry 


The geometry of a thin-walled, noncircular, 
cylindrical shell of length L and with an elliptical cutout 
located at the shell mid-length is shown in Fig. 1. The 
origin of the global Cartesian coordinate system, 
( x , y, z) is located at an end point of the longitudinal 


axis of the shell. As shown in Fig. 1, the x-axis 
coincides with the longitudinal axis of the shell. The y 
and z coordinates span the cross-sectional plane. A 
curvilinear coordinate system is also attached to the 
mid-surface of the cylindrical shell. The coordinates of 
points in the longitudinal, circumferential (tangential), 
and normal-to-the-surface (transverse) directions of the 
shell are denoted by (si, s 2 , ,S' 3 ), and the corresponding 
unit base vectors are {ei, e 2 , e 3 }. 

Following Romano and Kempner, 23 the non-circular 
cross-section of the cylindrical shell is defined as an 
oval with the coordinates, y and z expressed as 


y = R o X a J^ sin ~^- 

m= 1 , 3 , 5 , 7 ,... -^0 


and 



where £ represents the eccentricity of the oval cross 
section and R 0 is the equivalent radius of a circular 
cylindrical shell that has the same circumference as that 
of the oval cylindrical shell. The circumferential 
coordinate, s 2 varies between 0 and 27rR 0 . The 
derivation of Eqs. (1) along with the explicit forms of 
the coefficients a m (%) and b m (%) are given by 
Madenci and Barut. 25 

As derived by Romano and Kempner, 23 the 
coordinates y and z in Eqs. (la) and (lb) can be 
related to the radius of curvature of an oval-cross- 
section cylindrical shell, R(s 2 ,^) , by 


R{s 2 ,£) = 


l + £cos(2 s 2 /R 0 ) 


(2) 


Therefore, = 0 implies no eccentricity and 


corresponds to a circular cross section with radius R 0 . 
For positive values of the eccentricity parameter, £ , the 


z -coordinate becomes the major axis and the y - 


coordinate becomes the minor axis. For negative 
values of £ , the major and minor axes switch to the y- 


and z-axes , respectively. The range of values of the 
eccentricity parameter, , is bounded by — 1 < £ < 1 . 


As shown in Fig. 1, the cylindrical shell contains a 
cutout. The shape of this cutout is defined such that if 
the shell is cut along a generator and flattened into a 
plane, the cutout becomes an ellipse with major and 
minor axes denoted by a and b , respectively. For 
simplicity and convenience, the cutout is referred to 
herein as an "elliptical" cutout. Because the domain of 
the analysis shown in Fig. 2 corresponds to a similar 
flat region, a subsequent mapping of the ellipse to a unit 
circle is possible, which enables the use of Laurent 
series expansions for the local functions. Note that the 
special case of a "circular" cutout is given by a = b . 




Fig. 2 Computational domain of a cylindrical shell with an 
elliptical cutout. 
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In the flat analysis domain, the minor and major axes of 
the ellipse are aligned with a local coordinate system, 
(x,,x 2 ), whose origin is located at the center of the 
cutout and coincides with the origin of the parameter 
grid, given by constant values of s\ and s 2 . that forms 
the curvilinear coordinates (si, ^ 2 ) on the cylindrical 
shell mid-surface. The orientation of the elliptical 
cutout is arbitrary with respect to the longitudinal shell 
axis. Hence, the orientation of the local x,- axis (major 
axis) of the cutout and the longitudinal s, -axis of the 
cylindrical shell is denoted by the angle, iff . The 
elliptical coordinates, a and /? , representing a family 
of confocal ellipses and hyperbolas, respectively, are 
utilized in order to obtain the stress-resultant 
distribution in the direction tangent to the cutout 
boundary. The coordinate a is equal to 
a 0 = tanfr‘(Z>/a) on the particular ellipse that 
corresponds to the elliptical cutout. The other 
coordinate, , varying from 0 to 2 Jt , is known as the 

eccentric angle and is related to the (x, , x 2 ) coordinate 
system by x, = a cos (3 and x 2 = b sin /? . The eccentric 
angle /? is similar to the angle used for polar 
coordinates. 

The symmetrically laminated cylindrical shells 
considered herein are made of K specially orthotropic 
layers, and each layer has an orientation angle, 0 k , that 
is defined with respect to the s, -axis. Each layer also 
has elastic moduli E L and E r , shear modulus, G LT and 
Poisson’s ratio v LT , where the subscripts L and T 
represent the longitudinal (fiber) and transverse 
principal material directions, respectively. 

As for the shell thickness variation, the non-uniform 
wall thickness of the shell is denoted by h (s, , s 2 ) , and 
its variation is included by assuming that the thickness 
of each ply, t k , varies as a function of the curvilinear 
coordinates in the form 


f (^'n $2 ) f-n I 1 £ x Cos 


2 nm x s x 
L 


-£ 2 Cos 


where t k0 denotes the nominal thickness of the k" 
layer in the laminate, and the parameters ( m l , m 2 ) and 
( , £ 2 ) respectively, denote the wave numbers and the 


amplitudes of the periodic thickness variation in the 
longitudinal and circumferential directions. While the 
wall thickness of the shell is allowed to vary across the 
shell surface, the aspect ratio of the plies through the 
thickness is maintained, thus making the thickness 
variation of each ply to remain conformable to each 
other throughout the shell surface. A periodic thickness 
variation in the longitudinal direction is obtained by 


setting £ x =£ 0 and £ 2 = 0 , and in the circumferential 
direction by £ x = 0 and £-, =£ 0 . A shell with uniform 
thickness, t k = t k0 , is obtained by setting £ x = 0 
and £ 2 = 0 . 


Boundary Conditions and External Loads 

To facilitate a general imposition of prescribed 
boundary tractions, displacements, or rotations; the 
external as well as the internal edge boundary F of the 
shell is decomposed into 


r=r +r +r 

1 1 ( 1 ) ( 2 ) ( 3 ) 


As shown in Fig. 1, r (1) and r (2) denote the external 
edge boundary of the cylindrical shell and F ( 


represents the traction-free internal edge boundary 
around the cutout. The unit vector normal to an edge is 
represented by n . Throughout this paper, a variable 
with the superscript is treated as a known quantity, 
arising from the externally applied loads or from 
prescribed displacements and rotations. Also, the 
subscripts n , s , and t denote the directions normal, 
tangent, and transverse (through-the-thickness) to the 
boundary, respectively. The details of how prescribed 
edge loads and displacements are imposed in the 
analysis are presented subsequently. 


Prescribed edge loads 

External loads are applied to a shell by specifying 
values for the positive-valued stress resultants shown in 
Fig. 1. More precisely, the membrane loads applied to 
the ^ th boundary segment, T (<) , are given by 

Nu=< (5a) 

N n = t : (5b) 

where Nn and N 12 are the axial and shear stress 
resultants, respectively, defined in the cylindrical 
coordinate system. Likewise, shell-wall bending loads 
that are applied to the l* 1 boundary segment are given 
by 

M n = -»;* (6a) 

Mj u + 2 M l2 2 = t * - 2m* 2 (6b) 

where M u and M 12 are the pure-bending and twisting 
stress resultants, respectively, defined in the cylindrical 
coordinate system. Moreover, the left-hand side of Eq. 
(6b) is the Kirchhoff shear stress resultant of classical 
shell theory. 

As a matter of convenience, the analysis is 
formulated to also permit the specification of 
concentrated forces and moments that are transmitted to 
the ends of the shell as if through a rigid end-ring, as 
shown in Fig. 3. Presently, the concentrated force P n 
and the concentrated axial torque Ft are included in the 
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analysis. The force P* is simulated in the analysis by 
specifiying a uniform distribution of the axial 
displacement, with the unknown magnitude A n , such 
that 

J N n dT = Pl (7a) 

r if) 

Likewise, the torque F\ is simulated by specifiying a 
uniform distribution of the tangential displacement, 
with the unknown magnitude A s , such that 

J N n dr = P 5 * (7b) 

r (<> 

The analytical process that is used to ensure that the 
magnitudes of A n and A s correspond to the specified 

values of P“ and P * , respectively, is described in the 
following section and in Appendix A. 



Fig. 3 Application of a concentrated force through a rigid 
end-ring by using elastic springs along a shell edge 

Prescribed edge displacements and rotations 

Edge displacements and rotations are applied to a 
shell by specifying values for the displacements and 
rotations shown in Fig. 4 that correspond to the 
positive-valued stress resultants shown in Fig. 1. In 
particular, the axial and tangential displacements, 
k* and m* , respectively, that are applied to the £ A 
boundary segment, r (f| , are given by 
M,(n*e 1 )=M* (8a) 

M 2 [(e 3 xn)*e,] = M* (8b) 

Similarly, the transverse displacement m 3 and the 
rotation about an axis tangent to an edge P n that are 
applied to the I th boundary segment are defined by 
w 3 = u* (9a) 

m 3 , (n*ej ) = i?* (9b) 


As mentioned previously, these prescribed 
displacements are enforced through the use of elastic 
edge restraints (springs) to relax kinematic admissibility 
requirments on the functions that are used to represent 
the displacement fields. The uniformly distributed 
extensional and rotational springs that are attached to 
the shell edges in the normal, tangential, and transverse 
directions and used to enforce the kinematic boundary 
conditions are depicted in Fig. 4. 



Fig. 4 Types of unifonnly distributed elastic spring supports 
that can be prescribed along the edge of a cylindrical 
shell. 

Specifying appropriate stiffness values for the springs 
results in full or partial restraints along the shell edges. 
A zero value of the spring stiffness corresponds to a 
traction-free-edge condition. In contrast, a value of the 
spring stiffness that is large compared to the 
corresponding shell stiffness effectively corresponds to 
a prescribed zero-valued boundary displacement or 
rotation. This approach effectively yields a prescribed 
kinematic boundary condition in the limit as the relative 
stiffness of the spring becomes much greater than the 
corresponding shell stiffness. Similarly, values for the 
spring stiffness can be selected that correspond to a 
specified uniform elastic restraint along an edge, similar 
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to that provided by a rigid end-ring. This capability is 
important, and useful, because in some test fixtures or 
actual structures the edge supports may not be stiff 
enough to simulate a fully clamped boundary condition 
or flexible enough to simulate a simply supported 
boundary condition. 

As depicted in Fig. 4, the membrane displacements, 
u n and u s , and the transverse displacement, 

u t = « 3 along the t 4 boundary segment are restrained 
by extensional springs with stiffness values of S n , S , 
and 5, in the directions normal, tangent, and transverse 
to the boundary, respectively. In addition to the 
extensional springs, the edge rotations, rJ and , 

along the f th boundary segment are restrained by 
rotational springs with stiffness values of J n and J s 
that correspond to rotation about axes tangent and 
normal to the edge, respectively. 

Extensional springs in the directions normal and 
tangent to the shell edge, with stiffness values of s n 
and s s , are also used to simulate load introduction 
through a rigid end-ring, as shown in Fig. 3. Specifying 
values for the spring stiffnesses s n and s s that are 
relatively large compared to the corresponding shell 
stiffnesses causes the shell edge to behave as if a rigid 
end-ring is attached that produces the uniformly 
distributed displacements with the corresponding 
magnitudes A n and A s . The values for A n and A s 
that correspond to the specified concentrated loads are 
determined by using a penalty parameter approach. 
This approach enforces the difference between the edge 
displacements of the shell and the unknown uniform 
rigid end-ring displacements, (u n - A n ) and (u s - A s ) 
to vanish, while retaining the corresponding potential 
energy of the applied concentrated loads P* and P‘ . 

Kinematics and Stress-Strain Relations 

The kinematic equations used in the present study 
are based, to a large extent, on the assumptions of 
Love-Kirchhoff classical thin-shell theory. 
Specifically, the axial, circumferential (tangential), and 
normal (normal to the mid-surface) displacements of a 
generic point of the shell are denoted by U l (s l ,s 2 ,s 3 ) , 

U 2 (s j,s 2 ,s 3 ) and U 3 (s l ,s 2 ,s 3 ) , respectively. The 
corresponding displacements of a generic point of the 
shell mid-surface that share the same unit vector normal 
to the mid-surface are denoted by m, (.S', ,.s 2 ) , u 2 (s 1 ,s 2 ) 
and u 3 (.?, , s 2 ) , respectively. In classical shell theory, 
these displacements are related by 

(A j ^2 > ^3 ) — W 1 (Al j ^2 ) _ ^3 P\ (Al > ^2 ) ( 10 a ) 

U 2 (>$ 1 , s 2 , S 3 ) u 2 (.$) , s 2 ) s 3 fi 2 (^j , s 2 ) (1 Ob) 


C/ 3 (5i,5 2 ,53) = M 3 (5 1 ,5 2 )( 1 0c) 


where P 3 {s 3 ,s 2 ) and P 2 {s 3 ,s 2 ) are the mid-surface 
rotations about the s 2 and ,s’i axes, respectively, that are 
given by 


P\ (*$1 9 S 2 ) ^3,1 (*^1 9 S 2 ) 

(11a) 

P2 (^1 ’ *^2 ) — ^3,2 (^1 ’ *^2 ) W 2 CW 2 ) 

R(s 2 ) 

(lib) 


in which a subscript after a comma denotes partial 
differentiation. The corresponding linear membrane- 
strain-displacement relations are given by 


> 


u u 






( 1 ^ 

e 22 

> — < 

1 m 22 +-m 3 I 

Yu . 





( W l,2 “*“ W 2.1 J 


and the bending-strain-displacement relations are given 
by 


-u. 




K = { 


K~. 




R 


.2 J 


I ^ I u 3 12 u 2,1 


(12b) 


It is important to point out that the expression given for 
the change in surface twist due to deformation, K n , is 

that originally published by Love 26 ' 27 in 1888 for 
general shells, in terms of lines of principal-curvature 
coordinates, and derived in the book by Timoshenko 
and Woinowsky-Krieger 28 for circular cylindrical 
shells. As indicated by Bushnell, 29 the expression for 
K n vanishes for rigid-body motions in contrast to the 
corresponding expression presented in Reissner's 
version of Love's first-approximation shell theory (see 
Reissner, Kraus, andNaghdi ). Equations (12a) and 
(12b), and the more general forms presented by 
Bushnell, 29 are sometimes referred to as the Love- 
Timoshenko strain-displacement equations. 
Justification for this terminology is given by 
Chaudhuri. 33 

The stress-strain relations used in the present study 
are those of the classical theory of laminated plates and 
shells, 34 which are based on a linear through-the- 
thickness distribution of the strain fields. For a thin, 
symmetrically laminated cylindrical shell, with variable 
wall thickness, the relationship between the membrane 
and bending stress resultants and the membrane and 
bending strains is expressed conveniently in matrix 
notation by 
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N = A(Sj,s 2 )£ (13a) 

and 

M = D(jj,s 2 )K (13b) 

The membrane and bending stress resultants in Eqs. 
(13a) and (13b) are defined as 

K T ={N n ,N 22 ,N n } (14a) 

and 

M T ={M u ,M 22 ,M u ) (14b) 


It is important to reiterate that when shell-wall 
thickness variations are present, the membrane and 
bending stiffness matrices, A(5 1; 5 2 ) and D(5,,5 2 ),are 
dependent on the curvilinear surface coordinates s, and 


s 2 . 

It is convenient, in the present study, to combine the 
relations given in Eqs. (13a) and (13b) into the matrix 
form 

s = Ce (15) 

in which S, e and C are defined as follows: 


s r ={N r ,M r } 

e r ={c r ,K r } 


C = C ( Sl ,s 2 ) 


A(.s' |v s' 2 ) 0 

0 D(.s’ ]v s’ 2 ) 


(16a) 

(16b) 

(16c) 


Equations Governing the Response 

A general analytical approach for the exact solution 
of the equilibrium equations for a laminated-composite 
cylindrical shell with variable curvature is not 
mathematically tractable. Therefore, a semi-analytic 
variational approach that is based on the principle of 
stationary potential energy is used in the present study 
to obtain numerical results. Because elastic edge 
retraints are used as a means to relax the kinematic 
admissability conditions on the assumed displacement 
functions, and because a rigid-end-ring capability is 
used to impose shell-end force resultants, the potential 
energy consists of the elastic strain energy of the shell 
and the elastic edge restraints and the potential energy 
of the applied loads. In particular, the potential energy 
is expressed symbolically by 
*(q,A) = £/(q) + £l(q,A) + P(q,A) (17) 

in which U and il represent the strain energy of the 
laminate and the elastic edge supports (springs), and V 
represents the potential energy due to external boundary 
loads. Their explicit forms are presented in Appendix 
A. The symbol q is the vector of unknown, generalized 
displacement coefficients that arises from the 
mathematical representation of the mid-surface 
displacement fields that is used in the variational 
solution process. In particular, the mid-surface 


displacement fields are given symbolically by z/,(q), 
m 2 ( q), and « 3 (q) . The symbol A represents the vector 
of unknown edge displacements that arise from 
prescribing end loads. 

Subjected to the constraint equations that arise from 
the use of Lagrange multipliers, the equations 
governing the shell response are obtained by enforcing 
the requirement that the first variation of the total 
potential energy vanish. As discussed by McFarland et 
al., 35 because the constraint equations are not 
functionally dependent on spatial coordinates, s, and s 2 , 

the equations governing the response may be generated 
by modifying the total potential energy into the form 
n\q,A,l) = x(q,A) + W(q,l) (18) 

in which W is viewed as the potential energy arising 
from constraint reactions. In particular, 
tF(q,X) = I r Gq = 0 (19) 

where /. is the unknown vector of Lagrange multipliers 
and G is the known constraint coefficient matrix. 

Substituting the specific expressions for U (q) , 
U(q,A), F(q,A), and IF(qA) that arise from 

approximation of the surface-displacement field and 
enforcing the first variation of the modified form of the 
total potential energy to vanish lead to 

Sj f=Sq T [k qq q + S w q-s ?A A-f*-T* + G r k] 

+ SA t [s m A-S^ a q— P*] + J)/Gq = 0 (20) 

in which the matrix, k represents the stiffness matrix 

of the shell and requires evaluation of the 
corresponding integrand over a doubly connected 
region (see Appendix A for details). The spring- 
stiffness matrices, S and , are associated with the 

deformation of the shell edges and displacement of the 
rigid end-ring, respectively. The spring-stiffness 
matrix, s qA , captures the coupling between the 

displacement of the shell edges and the rigid end-ring. 
The vectors f*,T*, and P* arise from the prescribed 
boundary displacements, external tractions and 
moments, and the concentrated forces applied to a rigid 
end-ring, respectively. For the arbitrary variations 
(<5q, SA , andJk), the stationary condition requires 
that the following equations must be satisfied: 

[(k w +Sjq— S ?A A-f*-T* + G r k] = 0 (21a) 

[ s aa A-s^ a q-P* J = 0 (21b) 

G q = 0 (21c) 
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It is convenient, to express Eqs. (21a) - (21c) into 
the single matrix equation 

KQ=F (22) 


where K and F represent the overall, system stiffness 
matrix and the overall load vector, respectively. These 
matrices have the general, expanded form 


(23a, b) 



1 

K 

•C5 

•S5 

< 

<3- 

1/5 

1 

G r 


F*' 

K = 

T 

S AA 

0 r 

and F = < 

p* 


G 

0 

0 


0 


in which 

K <, =K +s „ 


and F =f* +T* 


(23c, d) 


The vector of unknowns, Q , that appears in Eq. (22) is 
defined as 


Q = 


q 

A 

X 


(24) 


Solving for the vector of unknowns in Eq. (22) yields 
all the information needed to obtain a complete 
variational solution to a specific problem. The accuracy 
of a solution depends on the number of terms included 
in the expressions for the local and global functions 
representing the displacement fields and converges to 
the corresponding exact solution as the number of terms 
increases. 


Displacement-field representation 

Representation of the mid-surface displacement 
field is a critical step in the variational solution to the 
problem. By relaxing the requirements for kinematic 
admissibility, the mid-surface displacement fields are 
represented in the present study by a combination of 
rigid-body modes, u Ri , and global and local functions, 

denoted by u i and u i , respectively; that is, 
u i =u Ri +u i +u i (25) 

where the values of the index are given by f = 1,2, and 
3. The rigid-body modes account for the overall or 
global translation and rotation of the shell, and are 
selected so that they produce neither membrane strain 
nor changes in shell curvature and twist. These terms 
are included for the completeness of the kinematics of 
the cylindrical shell. The presence of the appropriate 
displacement boundary conditions inherently eliminates 
the rigid-body motion. However, for cases where an 
insufficient number of kinematic boundary conditions 
are imposed, these rigid-body terms need to be 
eliminated, as discussed in detail in Appendix C. 
Following the complex-variable solution techniques 
used in the theory of elasticity, the local functions are 
expressed in terms of robust, uniformly convergent 
Laurent series (used for doubly connected regions) to 
enhance capturing steep stress gradients and 


deformations near the cutout. Complete sets of 
trigonometric expansions are used to primarily capture 
the overall global response of the shell. Here, 
completeness means that all the fundamental 
waveforms needed to construct the typical overall 
deformations of a shell are included in the set. 

For convenience, the displacement representations 
are rewritten in matrix form as 

u, = V 7 a R + V, r c,.+ f, 7 ' a (i = 1, 2) (26a) 

= Y r} a R + V 3 r c 3 + ff p (26b) 


An even more useful, compact form is given by 
u. =V. r q with f = 1,2,3 (27) 

where the vector of unknown displacement coefficients, 
q , is defined by 


q 7 ={a 7 ,c 7 ,c 7 ,c 3 ,a r ,p 7 } 


(28) 


In Eq. (28), the vector a R contains the unknown 
coefficients for the rigid-body motion of the shell, and 
the vectors a and p contain the real and imaginary 

parts of the unknown coefficients a, m and /?,„» , 
respectively, that are associated with the local 
functions. The vectors C, , where / = 1,2,3 , contain the 


real-valued unknown coefficients, c i( „ m) that are 
associated with the global functions. The explicit forms 
used herein for the unknown coefficient vectors a r , 


C, , a , and P that appear in Eqs. (26a) and (26b) 
along with the vector functions V, (and the 

corresponding subvectors V„, , V, , and V ; ) are given in 
Appendix B. 

In addition to the general representation of the shell 
surface-displacement fields, similar matrix expressions 
are needed for the displacements and rotations of points 
on the shell boundary. In the present study, the 
boundary displacement vector u r is introduced that 
consists of the mid-surface boundary displacements in 
the directions normal, tangent, and transverse to a shell 
edge, and the mid-surface rotations about axes that are 
normal and tangent to a shell edge. The boundary 
displacements in the directions normal, tangent, and 
transverse to a shell edge are denoted herein by w„ , u s , 
and u t , respectively. Similarly, the mid-surface 
rotations about axes that are tangent and normal to a 
shell edge are denoted by i? (I and if , respectively. In 
terms of the vector of unknowns defined by Eq. (28), 
the boundary displacements and rotations are expressed 
in matrix form by 

u r =Bq (29) 
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in which the boundary displacement vector, u r is 
defined by 

VL T r ={u n ,u s ,u„& n } (30) 

The matrix B is a known matrix of coefficients that is 
defined as 


B = 



(31) 


in which the sub-vectors, u T n , u j , uf and 0^ are 
known and defined by 


u* = (n*e,) V, (32a) 

= [(e 3 xn>e 2 ] V 2 (32b) 

uf = V 3 r (32c) 

and 

0f =(n.e 1 )V 3 r , (32d) 


Strain- and stress-resultant-field representation 

After defining the shell mid-surface displacement 
field in terms of the generalized coordinate q, the 
corresponding representation of the strains is obtained 
by substituting Eq. (27) into the strain-displacement 
relations given in vector form by Eqs. (12a) and (12b). 
This substitution yields 

£ = L f q (33a) 

and 


K = L v q 


(33b) 


where the strain-coefficient matrices L f and L s are 
defined as 

V r 

’ u 


L = 


L = 


y r + — V r 

y 2,2 ~ r y 3 

\ T +v r 

T 1,2 t 2,1 


~\ T 
” 3,1 1 

-V r + — V r + — v r 

“ 3,22 R 22 R 2 2 

-2V 1 + — V t 

^ T 3,12 1 n v 2,1 

K 


(34a) 


(34b) 


Next, the representations for £ and k are substituted 
into Eq. (15b) to obtain 

e = L q (35) 

where the overall strain-coefficient matrix L is defined 
as 

U=[U e I/J (36) 


Finally, the corresponding matrix representation of 
the stress resultants in terms of the generalized 
coordinates is obtained by substituting Eq. (35) into 
constitutive Eq. (15). The resulting vector of stress 
resultants is given by 

s = C Lq " (37) 


Constraint Equations 

In the generalized-coordinate representations for u, 
and u 2 , the coefficients c uo o> and c 2(0 o> associated with 
the global functions, iT, and u 2 , also correspond to 
rigid-body translation in the s', direction and rigid-body 
rotation about the s) axis, respectively. These two 
redundant rigid-body modes are eliminated by 
introducing constraint conditions using Lagrange 
multipliers. In particular, the unknown Lagrange 
multipliers A RRB{1) and 2 RSB{2) are associated with the 
redundant rigid-body modes. Also, multi-valuedness of 
the normal-direction displacement u 3 (s, , s 2 ) that arises 
from the presence of logarithmic terms in the Laurent- 
series-expansion for the local function — must be 
eliminated. The unknown Lagrange multipliers 2 sv(r) 

and 2 SV(s) are used herein to eliminate this multi- 
valuedness. Likewise, the rigid-body modes of the 
cylindrical shell must be eliminated by the Lagrange 
multipliers f RB( ;) ( j = 1,.., 6 ) if the specified kinematic 

boundary conditions are not sufficient enough to 
prevent them. In other words, the non-vanishing rigid 
body modes must be eliminated by introducing 
constraint conditions prior to the stress analysis in order 
for the overall system stiffness matrix K, given in Eq 
(22), to be nonsingular. 

These requirements on the representation of the shell 
displacement field are enforced by using constraint 
equations that use Lagrange multipliers. These 
constraint equations are functionally independent, 
forming a set of linearly independent equations equal in 
number to the total number of Lagrange multipliers. 
The Lagrange multipliers can be viewed as the 
reactions that are needed to enforce the corresponding 
constraints. In the present study, all of these constraint 
conditions are included in the matrix equation given in 
Eq. (19). The explicit form of the vector of unknown 
Lagrange multipliers, /. , and the known coefficient 
matrix, G , are given in Appendix C. 

Overview of Validation Studies 

A limited series of validation studies were conducted 
in the present study to determine the accuracy of results 
obtained by using analysis method presented herein. 
Specifically, the studies included circular and non- 
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circular cylindrical shells with either a circular or an 
elliptical cutout under uniform tension. The stress 
resultants around the circular and elliptical cutout for 
varying aspect ratios and orientations in a circular 
cylinder as well as the stress concentrations arising 
from a circular cutout in a non-circular cylindrical shell 
were computed. Comparisons of the stress-resultant 
distributions and magnitudes in the shells were made 
with the corresponding results obtained by using an in- 
house finite element program developed earlier by 
Madenci and Barut. 36 This finite element program has 
been validated, to a large extent, against previously 
published experimental and numerical results for stress, 
buckling, and post-buckling of thin-shell structures (see 
Madenci and Barut 37 ’ 38 ). Therefore, this finite element 
program is expected to serve as a reliable indicator of 
the accuracy of the analysis methods and results 
presented herein. Overall, the comparisons indicate 
very good agreement (less than 1% difference) between 
the corresponding results produced by the two analysis 
methods. For shells with high-aspect-ratio cutouts, 
differences of approximately 5% were obtained and 
found to be the result of insufficient mesh refinement in 
the finite element models. 

Selected Numerical Results 

Selected numerical results are presented in this 
section to demonstrate the utility of the analysis method 
presented herein and the potential for its use in 
developing design technology. These results elucidate 
the effects of loading condition, non-circular cross- 
section geometry, wall-thickness variation, cutout 
shape, cutout size, and cutout orientation on the 
intensity of stress-resultant concentrations near a 
cutout. Specifically, tension, torsion, and pure-bending 
loads are considered for 

[45° /- 45° / 90° / 0° / 90° /- 45° / 45° ] s quasi-isotropic 
shells with length L = 356 mm and made of graphite- 
epoxy plies. The nominal ply thickness 
ist k0 =0.14 mm, resulting in the total thickness of the 
shell given by h = 2.24 mm , and the ply orientation 

angles are measured with respect to the longitudinal 
shell axis. The Young’s moduli of each ply in the 
longitudinal, fiber direction and in the direction 
transverse to the fibers are specified as 
E l = 135. 0 GPa andis r = 13.0 GPa , respectively. The 
in-plane shear modulus and Poisson’s ratio of each ply 
are given by G LT = 6.4 GPa and v LT = 0.38 . 

The effects of varying the radius of curvature R 0 on 
the stress-resultant concentration along the contour of a 
circular cutout with radius a = 25.5 mm are shown in 
Fig. 5 for a circular cylindrical shell subjected to a 
uniform axial tension load. Four curves that correspond 


to values of R 0 /L = 0.5, 0.75. 1, and 1.25 are presented 
that show the tangential stress resultant, N # 

normalized by the far-field applied uniform stress 
resultant ,V 0 , as a function of position around the cutout 
(indicated by the "cutout angle", (Z >). As shown in Fig. 
5, the stress-resultant concentration is a maximum at 
<j)= 90° and 270° (at the net section of the shell) for 
each case and reduces from a maximum value of 
approximately 4.0 to a minimum value of 3.4 at the net 
section as the radius of curvature increases. In addition, 
the results show that the Y # (a,90°) stress-resultant 

concentration approaches the well-known value of three 
for an isotropic plate as the shell radius increases. Away 
from the net section, changes in the radius of curvature 
have a relatively small effect on the stress-resultant 
concentration. 



Fig 5 The effect of varying shell radius on the tangential 
stress resultants around a circular cutout in a quasi- 
isotropic circular cylindrical shell subjected to a 
uniform tension load. 

The effects of varying the circular-cutout radius on 
the stress-resultant concentration along the contour of a 
circular cutout is shown in Fig. 6 for a circular 
cylindrical shell with radius 7? 0 =381 mm and 
subjected to a uniform axial tension load. Five curves 
that correspond to values of the cutout radius a = 15, 
25.5. 30, 40, and 50 mm are presented that also show 
the tangential stress resultant N^{a,0 ) , normalized by 

the far-field applied uniform stress resultant, N 0 , as a 
function of the cutout angle <j) . The results in Fig. 6 
show that the stress-resultant concentration is a 
maximum at the net section of the shell for each case, 
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as expected, and changes significantly from a minimum 
value of approxiamtely 3.1 to a maximum value of 5.1 
at the net section as the cutout radius increases - an 
increase of approximately 65%. The results also show 
that the N^(a, 90°) stress-resultant concentration 

approaches the well-known value of three for an 
isotropic plate as the cutout radius decreases. Away 
from the net section, changes in the cutout radius have a 
much smaller effect on the stress-resultant 
concentration. 



Fig 6 The effect of varying circular-cutout radius on the 
tangential stress resultants around a circular cutout in 
a quasi-isotropic circular cylindrical shell subjected 
to a unifonn tension load. 

The effect of varying the elliptical-cutout aspect 
ratio, a/b , on the tangential stress-resultant distribution 
around the edge of a cutout in a cylindrical shell with 
radius R 0 =178 mm , and subjected to uniform tension 
is presented in Fig. 7. The orientation of the elliptical 
cutout is specified by i// = 0° . Two curves that 
correspond to the locations </> = 0° and 90° are 
presented that show the tangential stress resultant, 
N p j; (a n -P) normalized by the far-field applied uniform 

stress resultant ,V 0 , as a function of the cutout aspect 
ratio. As expected, the normalized stress-resultant 
concentration, (a 0 , ft) /N 0 , remains negative for all 

aspect ratios at </> = 0° , consistent with the expected 
Poisson effect, and the magnitudes are relatively 
insignificant at this location. In contrast, large stress- 
resultant concentrations are indicated at the net section 
(<j> = 90°) that diminish from a maximum value of 
approximately 17.0 for a widthwise, slot-like cutout 


with (a = 5 mm and b = 30 mm) or (a/b = 1/6) to a 
minimum value of 1.4 for a lengthwise, slot-like cutout 
(i a = 30 mm and b = 5 mm) or (a/b = 6) . 



Fig 7 The effect of varying elliptical-cutout aspect ratio on 
the tangential stress resultants around a cutout in a 
quasi-isotropic cylindrical shell subjected to a 
uniform tension load. 

The effects of varying the orientation of a high- 
aspect-ratio, slot-like elliptical cutout on the stress- 
resultant concentration along the cutout contour is 
shown in Fig. 8 for a circular cylindrical shell with 
radius 7? 0 =178 mm and subjected to a uniform axial 
tension load. The major and minor axes of the cutout 
are given by a = 30 mm and b = 5 mm , respectively. 
The orientation of the elliptical cutout, with respect to 
the longitudinal shell axis, is measured by the angle, iff . 
Three curves that correspond to values of if/ = 0°, 45°, 
and 90° are presented that show the tangential stress 
resultant at the cutout edge, N p /s normalized by the far- 

field applied uniform stress resultant, N 0 as a function 
of the cutout angle <j) . 

The results in Fig. 8 show that the stress-resultant 
concentration is the least pronounced for the case of 
iff = 0°. For this case, the cutout major axis is aligned 
lengthwise with the shell axis and the net section of the 
shell is the largest. The location on the cutout edge 
defined by (/) = 0° corresponds to where the edge of the 
cutout intersects the major axis. At this location, the 
edge of the cutout is in tangential compression 
(Npp/N 0 =-\.6), consistent with a Poisson effect. 

The location defined by <j> = 90° corresponds to where 
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the edge of the cutout intersects the minor axis; that is, 
at the net section of the shell. At this location, the edge 
of the cutout is in tangential tension (N^/N 0 =1.4). 

Between approximately ^ = 10° and 170° and between 
0 = 1 90" and 350°, the cutout width (and hence net 
section width) does not vary greatly. This attribute 
accounts for the corresponding flat regions in the 
iff = 0° curve shown in Fig. 8. 



Fig 8 The effect of varying elliptical-cutout orientation on 
the tangential stresses around the cutout in a quasi- 
isotropic cylindrical shell subjected to a uniform 
tension load. 

For the case of i// = 90° , the cutout major axis is 
perpendicular to the shell axis and the net section of the 
shell is the smallest. As before, the locations defined 
by 0 = 0° and 1 80° correspond to where the edge of the 
cutout intersects the major axis; that is, at the net 
section of the shell. The results in Fig. 9 show that the 
edge of this high-aspect-ratio cutout has extremely high 
stress-resultant concentrations at these locations 

( Npp /N 0 =17.) that have very step gradients. Between 

approximately 0 = 5 0 and 175° and between 

0 = 1 85° and 355°, the analysis predicts relatively 
benign variations in the stress-resultant concentration. 
The case of \j/ = 45° , exhibits stress-resultant 
concentrations that are, for the most part, bounded by 
the corresponding results for i// = 0° and 90°. The 
analysis also predicts very high stress-resultant 
concentrations where the cutout edge intersects the 
major principal cutout axis ( N^/N 0 = 8.2 ). 


The effects of varying the cross-section eccentricity 
(see Eq. (2)) of a tension-loaded oval shell with a 
circular cutout are shown in Fig. 9. The results in this 
figure correspond to the equivalent shell radius 
R 0 =381 mm and a circular-cutout radius given by 
a = 25.5 mm . Moreover, the tangential stress-resultant 
concentation at the shell net section, N^ia, 90°), 
normalized by the applied load N 0 , is shown as a 
function of the eccentricity parameter for the range of 
-0.15 <£ <0.15 . As indicated in the figure, negative 
and positive values of £ correspond to cylindrical 
shells with the largest cross-sectional width oriented 
parallel and perpendicular to the tangent plane that 
passes through the two points of the cutout edge that are 
on the surface generator that passes through the center 
of the cutout, respectively. A value of £ = 0 
corresponds to a circular cross-section and a value of 
£ = 0.15 corresponds to cross-sectional aspect ratio of 
0.9. 



Fig 9 The effect of shell cross-section eccentricity on the 
stress-resultant concentration in an oval 

quasiisotropic cylindrical shell with a circular cutout 
and subjected to a unifonn tension load. 

The results presented in Fig. 9 show that the stress- 
resultant concentration is affected benignly by the 
cross-sectional eccentricity. In particular, the stress- 
resultant concentration increases almost linearly with 
increases in the eccentricity parameter from 

A # (fl,90°)/A 0 = 3.5 to 3.6, which is slightly less than 

a 3% variation. This trend is understood by noting that 
the shells that correspond to negative values of £ are 
flatter near the cutout than those that correspond to 
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positive values of £ and, as indicated by the results in 
Fig. 5, are expected to have the lower values for the 
stress-resultant concentrations. 

The effects of longitudinal and circumferential 
periodic wall-thickness variations on the stress-resultant 
concentration at the net section of circular cylindrical 
shell with radius K,, =178 mm , circular cutout 
radius a = 25.5 mm , and subjected to uniform axial 
tension load are shown in Fig. 10. Two monotonically 
increasing curves that correspond to values of e x 
(with e, = 0 ) and e 2 (with £j = 0 ) are presented that 
show the tangential stress resultant N^(a, 90 °) , 
normalized by the far-field applied uniform stress 
resultant N 0 , as a function of thickness-variation 
amplitudes (see Eq.(3)) that range from 0 to 0.2. For 
the longitudinal thickness variation, the wave numbers 
used in Eq. (3) are m, =1 and m 2 = 0 . Similarly, for the 
circumferential thickness variation, the wave numbers 
used in Eq. (3) are m l = 0 and m 2 = 1 . 



Fig 10 The effects of longitudinal ( £j ± 0 and £ , = 0 ) and 
circumferential ( £j ± 0 and £^ = 0 ) wall thickness 
variations on the tangential stress-resultant 
concentration around a circular cutout in a quasi- 
isotropic circular cylindrical shell subjected to a 
unifomi tension load. 

The results shown in Fig. 10 indicate that the stress- 
resultant concentration at the shell net section increases 
as the magnitude of the thickness variation increases, 
for variations in either the longitudinal or 
circumferential direction. The maximum variation in 
the results is approximately 56%. Furthermore, the 
change in the stress-resultant concentration is slightly 


more pronounced for the circumferential thickness 
variation than for the longitudinal thickness variation. 
These increases are primarily due to a drastic loss of 
bending stiffness near the net section of the shell, as 
indicated by the wave numbers m x = 0 and nu = 1 , 
where the thickness of the shell near the center of the 
cutout is smaller. 

The effects of varying the radius of curvature R 0 on 
the stress-resultant concentration along the contour of a 
circular cutout with radius a = 25.5 mm is shown in 
Figs. 11 and 12 for a circular cylindrical shell subjected 
to a uniform torsion load and a pure-bending load, 
respectively. The pure-bending load corresponds to 
using t n = M 0 cos(/r/s 2 ) in Eq. (5a). Four curves that 
correspond to values of R 0 /L = 0.5, 0.75. 1, and 1.25 
are presented that show the normalized values of the 
tangential stress resultant as a function of position 

around the cutout. In Fig. 11, N go is normalized by the 
far-field applied uniform shear stress resultant, T 0 . In 
Fig. 12, Njf is nonnalized by the far-field applied 
uniform bending stress resultant, M 0 . 



Fig 11 The effect of varying shell radius on the tangential 
stresses around a circular cutout in a quasi-isotropic 
circular cylindrical shell subjected to a uniform 
torsion load. 

The results in Fig. 11 indicate that the stress- 
resultant concentration has identical maximum 
magnitudes at (j)= 45°, 135°, 225°, and 315° (at the net 
section of the shell) for each case, which corresponds to 
maximum diagonal tension and compression stress 
resultants associated with the shear stress resultants 
near the cutout. The magnitudes of the stress-resultant 
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concentration for these four locations reduces from a 
maximum value of 6.8 to a minimum value of 5.1 as the 
radius of curvature increases (33% variation). Away 
from these four locations, changes in the radius of 
curvature have a smaller effect on the stress-resultant 
concentration. The results in Fig. 12 indicate that the 
stress-resultant concentration for the shell subjected to 
the pure-bending load is quite similar to that presented 
in Fig. 5 for the corresponding tension-loaded shell. 



Fig 12 The effect of varying shell radius on the tangential 
stresses around a circular cutout in a quasi-isotropic 
circular cylindrical shell subjected to a pure- 
bending load. 



Fig 13 Stress resultant distribution near the cutout in a quasi- 
isotropic circular cylindrical shell subjected to a 
pure-bending load. 


Specifically, the stress-resultant concentration is a 
maximum at <j) = 90° and 270° (at the net section of the 
shell) for each case and reduces from a maximum value 
of 4.0 to a minimum value of 3.5 at the net section as 
the radius of curvature increases (14% variation). In 
addition, /M 0 approaches the well-known value of 

three for an isotropic plate as the shell radius increases, 
and away from the net section, changes in the radius of 
curvature have a relatively small effect on the stress- 
resultant concentration. For the case of R 0 /L= 0.5 
shown in Fig. 12, a contour plot of A # /M 0 near the 

cutout is shown in Fig. 13. The extent of the stress 
concentration at the shell net section ( (j) = 90° and 270°) 
is clearly captured by the analysis method presented 
herein. The highest stress-resultant concentration is 
N^(a,90°)/M 0 = 4 and it attenuates to the value of 

1.01 at a radius of about 80 mm (approximately three 
times the cutout radius), measured from the center of 
the cutout. 

Concluding Remarks 

A special-purpose, semi-analytical approach based 
on complex potential functions has been presented that 
can be used to investigate the behavior of thin, 
noncircular cross-section cylindrical shells made of 
laminated-composite materials and with a cutout, 
efficiently and parametrically. In particular, the effects 
of radius of curvature; elliptical cutout size, aspect 
ratio, and orientation; oval cross-section eccentrictity; 
wall-thickness variations; and loading conditions on the 
stress-resultant concentration near the cutout have been 
presented for a quasi-isotropic shell subjected to 
uniform tension, uniform torsion and pure bending. In 
addition, studies that were conducted to validate the 
analysis method have been described. 

A key finding of the results obtained with this 
analysis method is that the maximum tangential stress- 
resultant concentration near a circular cutout in a 
tension-loaded, circular, quasi-isotropic shell increases 
by approximately 18% as the shell radius-to-length 
ratio decreases from 1.25 to 0.5. Likewise, increases in 
the maximum tangential stress-resultant concentration 
as large as 65% have been found to occur with a five- 
fold increase in cutout radius. Results have also been 
presented that show extremely high tangential stress- 
resultant concentrations can occur for high aspect ratio 
elliptical cutouts whose principal axes are not aligned 
with the longitudinal axis of a tension-loaded shell. 

Additionally, results have been presented that show 
tension-loaded oval shells with a circular cutout on one 
of the flatter sides exhibit slightly lower tangential 
stress-resultant concentrations than the corresponding 
shell with the cutout on one of the more highly curved 
sides. Results have also been presented that show that 
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wall-thickness variations in either the longitudinal or 
circumferential directions significantly affect the stress 
concentration, with respect to that for the corresponding 
shell with a nominal thickness. The analysis also 
predicts that a quasi-isotropic shell with a circular 
cutout and subjected to pure bending that yields the 
maximum tensile stress resultant at the longitudinal axis 
of the cutout behaves similarly to the corresponding 
tension-loaded shell. The corresponding shell subjected 
to torsion was found to exhibit the maximum tangential 
stress-resultant concentrations at locations consistent 
with the maximum diagonal tension and compression 
near the cutout. Overall, the results demonstrate that the 
analysis approach is a powerful means for developing 
design criteria for laminated-composite shells. 
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Appendix A 
Strain Energy of shell 

Based on classical laminated shell theory, the strain 
energy of the shell can be expressed as 

U = — Js r e dA (38) 

in which A is the planform area of the shell mid- 
surface. Substituting the expressions for the resultant 
stress and strains, given in terms of the vector of 
unknown displacement coefficients, q, by Eqs. (35) 
and (37), leads to 

£/(q) = ~ j"q r (L r C L) q dA (39) 

A 

The matrix L involves the derivatives of the assumed, 
functional displacement representations, and C is the 
overall constitutive matrix defined by Eq. (16c). The 
expression for the strain energy is rewritten into the 
final form used herein as 


U{ q)=^q r k w q 

(40) 

where 


k „=(( L 'CL)^ 

(41) 


A 


The evaluation of this area integral is performed 
numerically by employing basic quadrature techniques. 
In this analysis, the quadrature points are pre- 
determined by employing standard triangulation of the 
entire domain as described by Shewchuk. 24 

Strain energy of elastic restraints 

The strain energy of the elastic edge restraints 
(springs), D , is expressed as 
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n = £ j s a (u a -u;) 2 dr 

1 t=l a=n,s,t r (0 

+ {ZZ J -lA^-^ fdr (42) 

Z t=l <*="•* r, () 

+ ^ZZ f ^(«a-Aj 2 ^ + 

2 <=i 

As depicted in Fig. 4, the boundary displacements 
u n ,u s , and u, along the I th boundary segment are 
restrained by extensional springs with the stiffness 
values S n , S, , and 5, , respectively. Likewise, the 
boundary rotations r? and d s are restrained by 
rotational springs with the stiffness values J„ and J s , 
respectively. 

In order to apply concentrated forces along the edge 
of a shell and introduce edge displacements that are 
similar to those introduced by a rigid end-ring or by the 
loading platens of a testing machine, additional springs 
are uses to simulate the load-introduction effects of a 
rigid end-ring. In particular, rigid-end-ring loads are 
introduced into the shell by using extensional springs in 
the directions normal and tangent to the boundary with 
corresponding stiffness values of and s s , as shown in 
Fig. 3. By specifying relatively large values for the 
spring stiffnesses s n and s s , the laminate edge behaves 
as if a rigid end-ring is attached that produces the 
uniform displacements A (l and A s . In contrast, a 
relatively small spring stiffness between the shell edge 
and the rigid end-ring eliminates the presence of a rigid 
end-ring. 

The desired form of the elastic-restraint strain 
energy is obtained in terms of the unknown vector q by 
substituting expressions for the boundary displacements 
and rotations, given collectively by Eq. (29), into Eq. 
(42). This step yields 

Z (q r si"q+n™;-2,'0 

Z i~\ a-n,s,t 

+ ^Z Z (q rj £q + iC)* -V 1 ^* ) (43) 

^ £-1 a=n,s 

\ 

s a Aldr-2q T ^A a 

where the matrices and represent the stiffness 
contribution of the extensional and rotational springs 
attached to the l th segment of the boundary. These 
matrices are defined as 

SL2= J S a u a u T a dr (, a = n,s,t ) (44a) 

r m 

and 



C= f ( a = n,s ) (44b) 

The matrix , representing the stiffness of the springs 
attached to the rigid end-ring, is defined as 

= f s u n a n T a dr ( a = n,s ) (45) 

4 

The load vectors, and , are associated with the 
prescribed boundary displacements and rotations and 
are defined as 

C* = J S a u* U„ dT ( a = n,s,t ) (46a) 


r « Y = J j a &a dT (a = n,s) (46b) 

r (0 

The vector, , is associated with the unknown end- 
displacements that correspond to a given concentrated 
load and is defined as 

= j x a u a dr ( a = n,s ) (47) 

The strain energies in the springs that arises from the 
known prescribed displacements ( u * , u* and tc * ) and 
rotations ( id* and id* ) are defined as 
&u(a) = J s a u a dY ( a = n,s,t ) (48a) 


n«> * = 

id(a) 


J J a $a dY ( a = n,s ) 


For convenience, the expression for the strain energy in 
the springs is recast in matrix form as 

Q(q,A)=iq r S w q + ^A r Saa A ^ 

q v a q r t a 

in which the matrices, S , and S qA represent the 

stiffness of the springs associated with the deformation 
of the laminate, the end-displacements and their 
coupling, respectively. These matrices are defined by 

= Z Z s «2 + Z Z + Z Z ( 50 a) 

£=1 a=n,s,t £=\ a=n,s £= 1 a=n,s 

>m = Diag[»;». *;•>, s ®,s;»]x2« 0 (sow 

v=[<" »?’ *'.' 1 »“] ( 50 c) 

The vector of unknown end-displacements, A , is 
defined by 

A r ={A« A?U?\A®} (51) 

The load vectors arising from all prescribed boundary 
displacements and rotations, f , is defined as 
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(58c) 


^ £=\ a=n,s,t 


,«) 


+\± z 

^ £-1 a=n,s 


(n 


(52) 


and the strain energy of all the springs due to prescribed 
displacements and rotations is 

^ £=\ a=n,s,t ^ £=\ a=n,s 


Potential of external loads 

The potential energy of the external tractions 
and f*) and moments (m* andw*) acting along 
the l * boundary segment, and the concentrated loads 
(P* and P* ) acting on the rigid end rings, is given in 
terms of the corresponding boundary displacements and 
rotations by 


r = -X Z f t' a u a dT 

e=i 

-z Z f K^dr-Z Z K K 


(54) 


t = 1 a=n,s r„ 


£-\ a=n,s 

Substituting the expressions for the boundary 

displacements and rotations, given in terms of the 
vector q, and combining terms in Eq. (55) yields 

F(q,A) = -q r r-A’p* (55) 

where the vector A , containing the uniform end- 
displacements A n and A s of the t lh boundary segment, 
is defined by 

A t ={A® A?\A« A?>} (56) 

The load vectors, T* and P* are defined by 

I = Z Z J C-I^r+XZ j ( 57a ) 

(57b) 


£-\ a=n,s,t r 


£=\ a-n,s y , 


and 


in which {l> P * , with ( a = n,s ) , represents the membrane 

forces applied on the boundary segment through a 
rigid end-ring. 


Appendix B 

Rigid-body modes 

As given by Madenci and Barut 24 , the rigid-body 
displacements ( u R] , u R2 and u R 3 ) of a cylindrical shell, 
defined with respect to the curvilinear coordinates, 
( Sj , s 2 , s 3 ), are expressed herein as 
u Rl =ot l -a 6 y + a 5 z 

u R2 = a 2 cos 0-a 2 sin 0-a 4 (ysinO + zcosO) 

+ a s x sin 0 + a b x cos 9 


(58a) 

(58b) 


u Ri = a 2 &md + a 2 cos 0 + a A (y cos <9 -z sin#) 

- a s x cos 8 + a 6 x sin 8 
where 8 denotes the angle between the radius of 
curvature at a point on the shell surface and z-axis as 
shown in Fig. 1. 


Global functions 

The global functions u i that are used to capture the 
overall deformations away from the cutout are 
expressed in terms of a series expansion of orthogonal 
functions of the form 

M m 

U, (z , ‘ s 2 ) = Z Z c,< m n) T m (z ) W„ ( S 2 ) (59) 

m= 0 w=0 

The symbols c i(mn) are the unknown real-valued 
coefficients, and T m (.s, ) and W n (s 2 ) are defined as 


TM) = 


C 


sin 


(m- 1 ) 


«■ + !) 


m = 0 
m = 1 

m > 1 


and 

W„(s 2 (8)) = 


(60a) 


(60b) 


cos(n#/2) n=0,2,4,6,8,-- 

sin((/? + l)#/2) n=l,3,5,7,9,-- 
in which -1 < £ < 1 and .S', is related to ^ as 
= ^T/2, with L being the length of the cylinder. 


Note that W n is periodical. These particular functions 

were chosen because they form a complete set of 
functions when used with Eq. (59). Flence, they are 
desirable for employing in energy based semi-analytic 
solution techniques such as the total potential energy 
principal that is used in this study. 


Local functions 

The local functions are expressed in terms of 
mapping functions that transform the contour of an 
elliptical cutout to a unit circle. These mapping 
functions are used permit the use of Laurent series 
expansions as local functions, which is desirable 
because Laurent series are analytic and uniformly 
convergent in domains with a circular hole. As a result, 
the use of mapping functions reduces the number of 
terms in the Laurent series significantly that are needed 
to adequately capture steep stress and strain gradients 
and local deformations near a cutout. In accordance 
with the principle of minimum potential energy, the 
local local functions are not required to satisfy the 
traction boundary conditions at the cutout boundary. 
Thus, the local functions, u . , are expressed in the form 
of Laurent series, in terms of complex functions, as 
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«j =2 Re 


u : = 2 Re 


u, =2 Re 


with 
P = \R 


±u" X 


w=l n——N 
n^O 


ix 2) z 


n= 1 n=-iV 

n*0 


"(/>) (61a) 


#(/» (61b) 


Y Y J3 F* (z ) 

/ j / j r' nm nm V Km ' 


m—\ n-—N 
n* 0 


77(A) (61c) 


2 + x 2 


(62) 

where the parameter N defines the extent of the 
complex series. In these series, a„ m and are the 
unknown complex coefficients that appear in Eqs. (26)- 
(28). The auxiliary function H{p) that defines the 
domain of influence of the local functions is expressed 
in a polynomial form as 
f f 


H{p) = 


1-10 


+ 15 


P_ 

\P„ 


-6 


r P_^ 

Pc 


0 < p < p 


0 (63a) 


l 0 P > Po 

with 

H(p o ) = H\p o ) = H\p o ) = 0 (63b) 

where the prime marks denotes differentiation with 
respect to the variable p and the parameter p„ denotes 
the radius of the region in which the local functions are 
effective. The purpose of chosing the auxiliary 
function is to prevent any possible linear dependency 
between the local and global functions and to restrict 
the influence of the local functions to a limited domain 
around the cutout. 

The complex functions up (z £m ) and m^ 2) (z em ) that 
appear in Eqs. (61a) and (61b) are defined as 


,0) 


( z cJ = cosi/f P m (z £m )~ sintp q m (z £m ) 


K ( z e„, ) = ^ Y P m ( z em ) + cos yf q m (z £m ) 


(64a) 

(64b) 


where the complex constants p m and q m are given by 
Pm = «1 1 PL + a n - PePem (65a) 

4m = a n Pern + «22 ' Pen “ «26 (65b) 

In Eqs. (65a) and (65b), the unknown complex 
constants, p rm , are the roots to the characteristic 
equation associated with membrane deformation, i.e., 

(66) 


a \ \Pem 2 a \bPem + (^26 + °66 )Pem 


^ a 26 Pern + a 22 b 

in which the coefficients a tJ are the coefficients of the 

flexibility matrix a , which is the inverse of the stiffness 
matrix A defined by Eq. (13a). Both the flexibility and 


the stiffness matrices, a and A , are measured with 
respect to the local coordinate system (jq ,x 2 ) . The 
angle, \j/ represents the orientation of the local 
coordinate system with respect to the global coordinate 
system, (.s', , .s 2 ) . 

The complex potential function, (z £m ) , 

appearing in Eqs. (61a) and (61b) is defined as 

®L( Z em) = ZL ( 67 ) 

in which the mapping functions, % £m , map a cutout onto 
a unit circle. The mapping functions for an elliptical 
cutout, introduced by Lekhnitskii 22 , are given by 

_i_ F 2 2 2 7,2 

Z ~, ± \l Z em- a -Psm b 


Zem = ~ 


(m = 1, 2) 


(68) 


0 -iPlnP 

where z £m = x 1 + p £m x 1 , a and b are the major and 
minor axes of the elliptical cutout, and i = \f— 1 . The 
sign of the square-root term is chosen so that \p £m | > 1 

(i.e., the mapped point is guaranteed to be on or outside 
the unit circle). 

Inverting the mapping function provides (O em (p £m ) 
as 


z =o) (£ l = r E - — 

Em fm \7frn / em^em r 


L 


(69) 


in which 

r em=\{a -iP em b ) , s em =^(a +ip £m b ) (70a, b) 

The unknown complex constants p eX and ju e2 , and their 
complex conjugates, i.e., p ei = p r[ and ,U f4 = p £l , are 
the roots obtained from the characteristic equation 
associated with membrane deformation 

The complex potential functions, F* nm (z Km ) in E. 
(61c) are defined as 


Pfim ( Z k'm ) 


pn ^ Kin pn-2 . i 

+ , - 'Km ’ r 2 r ^ 

n + 1 n — \ 


-fZL-Srm ln^„ n = 1 
«=- 1 


(71) 


in which the expressions for the mapping function t; Km 
and the constants r Km and s Km have the same form as the 
corresponding expressions for p £m , r em , and s em given by 
Eqs. (68) - (70a, b), except that the subscript £ is 
replaced by k . 

The complex variables z Km are defined by 


Z Km T PenPl 


(72) 
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in which the unknown complex constants fl Kl and li k1 
and their conjugates, i.e., ll k , = Jl k l and /t.. 4 = Ji K2 , are 
the roots obtained from the characteristic equation 
associated with the bending equilibrium equation 
D 22 juL + 4 D 2 X m +(20 12 + 4 D 66 )/4 (?3 

+ 4 D l6 n rm +D n = 0 

where Z> are the components of the bending stiffness 

matrix D (see Eq. (13b), which is defined with respect 
to the local coordinate system, (x, ,x 2 ) . 

It is important to note that the local functions in Eq. 
(61) satisfy the in-plane and bending equilibrium 
equations of a homogeneous, flat laminate of uniform 
thickness, not a cylindrical shell. Therefore, the roots 
to the characteristic equations, Eqs. (66) and (73) serve 
as approximation to their exact values which are not 
mathematically tractable. Because the solution 
procedure is based on the principle of minimum 
potential, their exact values are not necessarily 
required. However, they capture the stress 
concentration and local deformation near the cutout in 
cylindrical shells because these functions possess the 
inherent solution characteristics. They satify the 
equilibrium equations exactly as the radius of curvature 
approaches inifinity and are uniformly convergent in a 
doubly connected region. 

In the displacement representations defined by Eqs. 
(26a) and (26b), the vectors, a R , c, , a , and P are 


defined as 

a R ~ {^Rl> a R2^R3’ a R4^ ^R5’ ® Rb\ (74a) 

c ; = { c i(oo)> c i(io)> c i(oi)>—> c i(MO) (74b) 

’ C i((M-l)l)’”"’ C i(l(M-l))’ C <(OM)J 

« r ={«-« ,a^ +1 ,....,a(,,af, } (74c) 

in which 

<={< 1 X 2 } (75) 

with 

a '„j ={ R e [«„,], Im [«„,-]} (76a) 

and 

p r ={p^ ,p^ +1 „. ..,p! 1 ,p^„...,p;_ 1 ,p 5, } (76b) 

in which 

Pi ={PI„PI 2 } ( 77 ) 

with 

Pl y ={Re[^], Im [ft]} (78) 


The vector functions, V, , with i = 1, 2, 3 associated 
with the unknown generalized coordinates q appearing 
in Eq. (27) are defined as 


v;={v' l .v’,o r .o',v l '.o') 


'M 

'M' 


V J 0' V r O r V 


VI, ,o T ,o r ,v, r ,o r 


, r .5 r } 

•V/} 


(79a) 

(79b) 

(79c) 


where 


O r = {0,0,0,.., 0} of order [jM+l)(M+2)/2j (80a) 

0 7 ={0,0,..,0} of order 8 N (80b) 

The vectors associated with rigid body motion are 
Vjj ={l,0,0,0,z,-y} (81a) 


\ ds 2 ds 2 


dz dy 

y z - 

ds, ds 1 


dz dy 

,—x ,—x — 

ds, ds 2 


(81b) 


VJ = 0 .-— 

\ ds 2 ds 2 


dy dz 

+ z — 
ds, ds 




2 j 


(81c) 


dy dz 

,—x ,—x 

ds , ds 2 

Similarly, the vectors associated with the global 
functions are 

V = {T 0 (sJW 0 (s i ms l )W 0 (s 2 \ 

To ($1 )Hj ( s 2 )> T 2 ( 5 , )W g (s 2 ), 

T^W^XT^WM.. (82) 

...JMW^XT^w.is,),.. 

..., T x (s t )W M _ l (s 2 X T 0 (s 1 )W M (s 2 )} 
in which the expressions for T’(sj) and IV : (s 2 ) are 

given by Eq. (60), and 
V T ={v r V T V T 

= , ( 83 ) 

with 

V T =\\ T X T 1 (841 

y i(n) ( y i(«l)’ y ;(« 2 )J 

in which 

% nj) ={2Re[tt®4»;] -2Im[ M ®<I>;]} (85) 

where i, / = 1,2 and 


V 7 =(v T V r V r 

T 3 ( r 3(-W) > t 3(-W +!)>••••! t 3(-1) 


y T y T 

> ” 3d) 5 * **’ ’ ; 


'3(1)’ 

with 

y T = {y r y T 1 

y 3 (h) ( y 3(»l)’ y 3(«2)J 

in which 


V 

’ y 3( 


3 (N -1)’ ’ 3(Af ) 


( 86 ) 


( 87 ) 
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V 3 w) = {2 Re [F n * ] , -2 Im [F b * ] J , 0 = 1,2) (88) 


Appendix C 
Constraint Equations 

The unknown vector of Lagrange multipliers, X and 
the known coefficient matrix, G , in Eq. (19) are 
defined by 

^ — { / ^RRB(1)> / ^'RRB(2)> AsF(r), / ^ir(i), / ^iS( 1 ),-", / ^«B(6)} (89a) 


Or 

§RRS(1) 

o r 

¥ 

If 

<F 

Or 

o r 

§RRS( 2) 

o r 



Or 

o r 

¥ 

¥ 

if 

g SV(s> 

Or 

Or 

¥ 

o r 

if 

g W) 
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o r 

¥ 

o r 

if 

¥ 

§RB(2) 

¥ 

¥ 

o r 

if 

¥ 

§RS( 3) 

¥ 

¥ 

o r 

if 

¥ 

§RR( 4) 

¥ 

¥ 
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¥ 

¥ 

§RS( 5) 

¥ 

¥ 

¥ 

if 

¥ 

§RR( 6) 

¥ 

¥ 

¥ 

if 

¥ 


(89b) 


in which the constant coefficient vectors, 
Srksii) an d§iu(B( 2 ) are associated with redundant rigid- 
body modes, g sr(r) and g srw with single valuedness of 


the radial displacement component, and g /(;j(/1 with the 

rigid-body modes introduced by the global functions 
defined in Eq. (59). The rigid-body modes must be 
eliminated in the absence of a sufficient number of 
specified kinematic boundary conditions. These terms, 
as well as, the vectors with zeros are defined in the 
following subsections. 


Redundant rigid-body modes 

The coefficients, c 1(00) and c 2(00) in Eq. (59) for the 

global displacement functions produce additional rigid- 
body translation in the S) direction and rigid body 
rotation about the .v axis, respectively. Because these 
rigid-body modes are already represented by a R and 
a R6 in Eq. (58), the redundant rigid-body motion 
arising from the presence ofc 1(00) and c 2(00) , must be 

eliminated in order to obtain a unique representation of 
the displacements. These redundant rigid-body modes 
are eliminated by using the constraint conditions 

2rrB(1) C 1(00) = 0 (90a) 

^RRB( 2) C 2(00) — 0 (90b) 


in which the unknown Lagrange multipliers are denoted 

by 2 /(/(;;(l) and 2 /(/(/j(2) . In terms of the vector of 

unknowns, q , these constraints are rewritten in vector 


form as 

V 

Cl 



a 


IP 


(91a) 


(91b) 


in which the vectors of zeros, 0^ , 0 2 and 0 2 are 
defined as 

0 2 = {0,0, 0,0, 0,0} (92a) 

0 r ={0,0,0,. .,0} of order [(M+l)(M+2)/2] (92b) 

0 r ={0,0,0,. .,0} of order 8N (92c) 


The constant coefficient vector, g RRB is defined as 
g T RRB = {1,0,0,. .,0} of order [(M+l)(M+2)/2] (93) 


Single-valuedness of the radial displacement 
component 

The multi-valuedness of the normal displacement 
component that arises from the presence of logarithmic 
terms in the local expression for the radial displacement 
component in Eq. (61c) must be rendered single valued 
in order to obtain a unique solution. The logarithmic 
terms associated with complex constants, J3_ lm and j3 lm , 
with(»; =1,2), in Eq. (71) result in two real constants. 
Representing the complex variable of the Laurent 
series, <f m = p m e' e "‘ , the single-valuedness requirement 
is enforced as 

u, (L = pj 6 ■ ) - u z (£, = pj ^ +2 *’ ) = 0 (94) 

Associated with the complex constants, 
f_ Xm and P hn with (m = 1, 2) , in Eq. (71), this condition 
yields 

2R e{i{K n (L=P m e W ") 

(95a) 

-K m (Z m =pJ ie - +2 ’ > )}0-i m } = o 

or 
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Im Zk A i,„] = () 


and 


2Re{t{k m (4 =Pj e ") 

m = 1 

-^;(^=Ak (<, ” + ^)}Aj=o 


or 


Im Zt^ A„]=o 


(95b) 


(96a) 


(96b) 


In order to ensure single-valuedness, these constraints 
are enforced as 

2 


4 skm Im 


and 


4 sk« Im 


Zk™ A-,„ 


Zk™ A, 


= 0 


= o 


(97a) 


(97b) 


in which the unknown Lagrange multipliers are denoted 
by 

Av(r) and Av(s) . 

In terms of the vector of unknowns, q, these 
constraint conditions can be recast in matrix form as 

fa. 


V, {oLo r ,o r ,o r ,o r ,g^ (r) } 


Arc,, {o«,o r ,o r ,o r ,o r , gk (s ) } 


c, 

«2 

C 3 

a 

P 

a R 

c, 

c 2 

« 3 

a 

Pj 


= 0 


(98a) 


= 0 


(98b) 


where the constant coefficient vectors g SF(r) and g SF(s) 
are given by 


JT _ f T T T 

~ \8sr (/)(-#) ’ 05F(yX-iV+I) 8. 


*SV(j) 


SV(j)(- 1) 


T T T 1 

’ 8sr( /xi) 8sr( i)(N- 1) ’ 8 sv( i)(N) J 


SV(j)i 1) ’••••’ &SVUXN-1) ’ &SV(j)(N) 

in which 

gJrirxi) ={ Im k],Re[r 1 ], I m[r 2 ],Re[r 2 ]} 

= {0,0, 0,0} if «=£l 


&SV(r)(n) 

and 


Skos-ifi, = {lm[s, ] , Re [s, ] , Im[s 2 ] , Re [s 2 ]} 
= {0, 0, 0, 0} if /7 ^ 1 


gsF(j)(») 

with -N <n<N . 


(99) 

(100a) 

(100b) 

(100c) 

(lOOd) 


Rigid-body modes 

In the absence of kinematic boundary conditions, 
the rigid-body modes of the displacement field are 
eliminated by enforcing the constraint conditions in the 
form 

^■RiHjAmj) = 9 (101) 

where j = 1, 2, 3, 4, 5, 6 , and the unknown Lagrange 
multipliers are denoted by A s/J(/| . In terms of the vector 
of unknowns, q , these constraints are rewritten in 
vector form as 

fa. 


^ O ){gLo-,,0 r ,0 r ,0 r ,0 r ,0 r } 


= 0 


(102) 


are 


where the constant coefficient vectors, g RB(J) 
defined as 

gL 0 )=k AjAjAjAjAj) (io3) 

in which 8 is the Kronecker delta. 
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